An Embedding Technique for the Solution of Reaction-Diffusion 
Equations on Algebraic Surfaces with Isolated Singularities 

Parousia Rockstroh • Thomas Marz • Steven J. 
Ruuth 



Manuscript as of November 26, 2012 



Abstract In this paper we construct a parametrization-free embedding technique for numerically 
evolving reaction-diffusion PDEs defined on algebraic curves that possess an isolated singularity. In 
our approach, we first desingularize the curve by appealing to techniques from algebraic geometry. 
We create a family of smooth curves in higher dimensional space that correspond to the original 
curve by projection. Following this, we pose the analogous reaction-diffusion PDE on each member 
of this family and show that the solutions (their projection onto the original domain) approximate 
the solution of the original problem. Finally, we compute these approximants numerically by 
applying the Closest Point Method which is an embedding technique for solving PDEs on smooth 
surfaces of arbitrary dimension or codimension, and is thus suitable for our situation. In addition, 
we discuss the potential to generalize the techniques presented for higher-dimensional surfaces 
with multiple singularities. 

Keywords Closest Point Method • implicit surfaces • surface-intrinsic differential operators • 
Laplace- Beltrami operator • blow-up • singularity • resolution of singularity • singular differential 
operator 



1 Introduction 

In this paper we study the problem of evolving a reaction-diffusion PDE on algebraic surfaces that 
contain isolated singularities. Singular problems have been studied extensively from an analytic 
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perspective within the differential geometry community, originating with the work of Cheeger in 
[5], where a functional calculus was given for the Laplace operator on a cone. This in turn led to 
fundamental solutions of the heat equation and wave equation on a cone as given in [G] and [7]. 
More recently, researchers such as Jeffres and Loya in [18] have derived asymptotic expansions for 
solutions of geometric PDEs on conical manifolds. Our goal is to construct a parametrization-free 
embedding-based numerical method that approximates geometric PDEs posed on surfaces with 
cusp singularities. In addition to being an interesting problem in its own right, the construction 
of such a numerical method might be applied to the approximation of PDEs on smooth but 
highly curved surfaces such as the simulation of heat flow on filaments of high curvature or the 
development of cortical maps in which regions of high curvature occur along the sulcal lines. 

Our approach will be to use and extend existing embedding techniques. Within the past decade 
several numerical embedding methods have been created for solving PDEs that are posed on 
smooth surfaces. These include level-set methods such as those presented in [15] and [14] as well 
as the Closest Point Method as presented in [23] and [2G] and further extended in [20] and [21]. 
However, the existing embedding techniques are low-order accurate or inconsistent when applied to 
PDEs on surfaces with singularities. This is a result of the fact that within the class of embedding 
techniques, a smooth representation of the surface is needed in a narrow computational band 
surrounding the given curve or surface. When the curvature of the surface increases it becomes 
more difficult, and computationally expensive, to discretize in such a way that the region of high 
curvature is accurately captured. Moreover, when the curvature becomes infinite, as is the case 
with a cusp, classical numerical embedding methods may become inaccurate or fail due to the loss 
of smoothness in the surface representation. 

Our approach to evolving PDEs on singular surfaces is to evolve a modified PDE on a reg- 
ularized version of the surface. The first step here is to resolve the singularity in the underlying 
domain. To this end we employ a standard procedure from algebraic geometry which is known 
as "blowing-up" [17,16,27]. We construct the blow-up map in such a way that it produces a 
one-parameter family of smooth surfaces that approximate the original singular surface. The reg- 
ularization leads to the same reaction-diffusion PDE on each of the smooth surfaces within the 
one-parameter family. Solving these PDE problems yields accordingly a family of functions which 
converges to the solution of the original problem posed on the singular surface. 

By this approach, the domain of the PDE will change with the parameter. This is an effect 
which is disadvantageous for the numerics. Thus, we transform the PDE problems so that they 
all have the same smooth surface as a domain. This produces a one-parameter family of variable 
coefficient PDEs with coefficients depending now on the parameter which was introduced by 
the desingularization procedure. The resulting equations are now ready for the application of an 
embedding technique using standard uniform grid numerical methods on a smooth domain. 

In the construction of our numerical technique, we must be mindful of the fact that the blow-up 
procedure generally embeds the regularized surface into a higher dimensional space. It is therefore 
necessary to choose a numerical method that is effective for arbitrary co-dimensional embeddings. 
The Closest Point Method is one such embedding method, as shown in [21] and [26], and is our 
method of choice. 

The paper is organized as follows. In Section 2, we review the definitions relating to smooth 
embedded surfaces and surface-intrinsic differential operators. Section 3 introduces our problem of 
interest: a simple reaction-diffusion equation posed on the closed planar curve given by y 2 — x 3 —x 4: . 
This curve has a single cusp singularity at the origin. In addition, wc construct an analytical 
solution to the given problem which will be used later to assess our numerical results. In Section 4 
we describe the regularization of the problem in detail. In particular, we demonstrate the blow- 
up procedure, construct analytical solutions to the regularized problems, prove the convergence 
to the solution of the original reaction-diffusion equation, and describe the transformation which 
prepares the problem for an embedding technique. In Section 5, we review the basics of the Closest 
Point Method and apply it to the regularized problem. Moreover, we give convergence studies that 
demonstrate the robustness of the method. Section 6 discusses the extension of the presented 
method to general algebraic surfaces with finitely many singularities. 
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2 Algebraic Surfaces and Surface Intrinsic Differentials 

This section begins by collecting the definitions relating to smooth embedded surfaces and surface- 
intrinsic differential operators. 

In this paper, we consider real algebraic manifolds, which is precisely the class of curves and 
surfaces S C R" given implicitly as the solution of a system of m = codim S polynomial equations 
over the field R. We write this system as 

cp(x) =0 « xeS (1) 

with <p : R ra — > R m , where each component ipi, i G {1, . . . , m}, is a polynomial in n variables Xj, 

3 G {!,•■■,«}• 

Definition 1 (Tangent space) The straight line x + iv, t € R is tangent to the surface S at x e S 
if the root at t = of the polynomial y(x + iv) (a polynomial in t) has a multiplicity of at least 
2. By Taylor expansion 

( / 5(x + fv)=tD(^(x)-v + t 2 p(t), (2) 

where t 2 p(t) is the polynomial remainder and p is some polynomial, we see that this multiplicity 
requirement is equivalent to 

LV(x) ■ v = . (3) 

In this case we call v a tangent vector of S at x, and v is in the tangent space T X S. Moreover, by 
(3) the rows of D(p(x) span the normal space AT X S. 

Note that in the first part of Definition 1, we define tangency without the use of limits or deriva- 
tives, as is standard practice in algebraic geometry, see e.g. [27]. This will aid us in characterizing 
and resolving singularities as we will see later in the paper. 

Definition 2 (Regularity/Singularity/Smoothness) Let S C R™ be a real algebraic surface, then: 

(a) a point x £ S is called regular if D(p(x), the Jacobian at x, is a full rank matrix; in this case 
the columns of Dcp(x) form a basis of the normal space Af x S, 

(b) a point x g S is called singular if Dtp(x) is rank-deficient, 

(c) the surface S is called smooth if all points x g S are regular. 

Given a regular point x e S we set N(x) := D(p(x) T , hence N(x) is the matrix that contains the 
normal vectors as given by the implicit description. Moreover, we denote the orthogonal projector 
that projects onto the tangent space T^S of x € S by P(x). For regular points x£S one can write 
the matrix P(x) as 

P:S^R nX7 \ P(x) =/- N(x) • N(x) t , N(x) = £>y(x) T , (4) 

where N(x)t := (N(x) T N(x)) _1 N(x) T denotes the pseudo-inverse of N(x). In this paper we will 
be particularly interested in cases in which surfaces posses isolated singular points, hence N is 
full-rank and P is given by (4) almost everywhere. 

We now proceed to define a calculus for surface functions that does not depend upon parametriza- 
tions. Such a calculus has been introduced earlier (see for example [1,8,10,13,21,24]); here we will 
use much of the same notation and terminology as [21]. We point out here that throughout the 
paper surface functions, i.e., functions of the form / : S —> M m , are not required to be polynomials. 

Definition 3 Let S be a real algebraic surface and x G S be a regular point. We consider the 
following C^-smooth surface functions: a scalar function u : S — > R, a vector-valued function 
/ : S —> R m , and a vector field g : S —> R™ with C 1 -smooth local extensions into R™ called ue, 
!ei 9e (see [21]). We define the surface gradient Vs, the surface Jacobian D$, and the surface 
divergence divs at the regular point x by: 

V 5 u(x) T := V W£ (x) T ■ P(x) , 

D s f(x) := Df E (x) ■ P(x) , 
div S 5(x) := trace(D s g(x)) = trace(D# B (x) • P(x)) . 
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Here V is the gradient and D the Jacobian in the embedding space K™ applied to the extensions 
of the surface functions. 

By combining the above operators, we may form higher order differential operators on surfaces. 
One such operator is the Laplace-Beltrami operator, i.e., the Laplace operator for surfaces, which 
we define below. 

Definition 4 Given a smooth surface 5, the Laplace-Beltrami operator is defined by: 

A s u(x) := div s (V5u)(x), (5) 
for C'-smooth, I > 2, scalar surface functions u. 

Remark 1 In some literature (e.g. [3], [4], [19], [11]) the Laplace-Beltrami operator is defined by 
Asu = — divs(Vsit) to give a positive operator. This paper defines the Laplace-Beltrami operator 
to be a negative operator, as is standard when studying diffusion processes on surfaces [25]. 



3 Reaction-Diffusion on a Cuspidal Curve 

The differential operators given in Section 2 are well-defined on smooth regions of a surface. We 
are, however, interested in examining the case where the surface S posses an isolated singularity. 

For now, we focus on a sample problem which consists of a simple reaction-diffusion equation 
posed on a closed curve that has an isolated cusp singularity at the origin (x,y) = (0,0), but is 
smooth everywhere else. This section also constructs a parametric solution to the sample problem 
which we use later in the paper to validate our results. 



3.1 The Reference Problem 

The domain of our reference problem is the algebraic curve given by 

y 2 = x 3 - x 4 . (6) 

Formally, we let ip : M 2 — s- R be given by ip(x, y) — y 2 — x 3 + x A so that 

S = {(x,y)eR 2 :^(x,y)^0}. (7) 

For the remainder of the paper, unless stated otherwise, we will take this as our definition of S. 
Figure 1 shows a plot of S. Clearly, the origin is a singular point. We verify this using Definition 2: 
the Jacobian of if is given by Dip(x, y) = (— 3x 2 + 4x 3 , 2y). Since this is a 1 x 2-matrix, the singular 
points are those where the Jacobian vanishes, D<p(x,y) — 0. We find that the only singular point 
is the origin (x, y) = (0, 0). Consequently, the set S \ {0} is a smooth curve by Definition 2(c). 
We now state the reference problem. Consider the reaction-diffusion problem 

d t u = A s u - pi 2 u , t g (0,oo) , (x,y) € S\ {0} (8a) 

posed on the smooth part of S, i.e., on S \ {0}. By excluding the origin we introduce an interior 
boundary. The corresponding "gap" is closed by imposing the following boundary conditions 

u(t,0+) =u(t,0_) , (8b) 
d T u{t,0 + ) = <9 T u(i,0_) , ie(0,oo). (8c) 

Condition (8b) is a continuity requirement while (8c) is Kirchhoff 's circuit law regarding the fluxes 
into the interior boundary node. Here, T is the counter-clockwise tangent field of the curve 5\ {0} 
and dru(t,x,y) = T T (x,y) ■ V su(t, x,y) is the corresponding directional derivative. The positive 
and negative subscripts denote one-sided limits. The initial condition is given by: 

u(0,x,y) = u (x,y) , (x,y)eS\{0}. (8d) 

The initial boundary value problem (IBVP) specified by equations (8a)-(8d) gives our reference 
problem. 
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Fig. 1: The red curve is given by y 2 = x 3 — x 4 and possesses an isolated singularity at the origin. 



3.2 Parametric Solution of the Reference Problem 

We will now construct a parametric solution to the reference problem. This exact solution will be 
used to perform numerical convergence studies on our method. 

Let 7 : (0, L) — > S \ {0}, s — > j(s), denote a counter-clockwise arc-length parametrization 
of S \ {0} with lim 7(s) = = lim 7(s). The counter-clockwise tangent field is parametrized 

s— >0 s— >L 

by T o 7(s) = t'(s), and one-sided limits of a surface function / are given, in terms of 7, by 
/(0_) = hm/o 7 ( s ),/(0 + )= lim/o 7 ( s ). 

s— >0 s— s-L 

We denote by u(t, s) := u(t, 7(s)) the parametric solution of the reference problem correspond- 
ing to the arc-length parametrization 7. Then, u must satisfy the PDE: 

d t u = d 2 s u - n 2 u , t € (0, 00) , s G (0, L) (9a) 

with periodic boundary conditions 

lim u(t, s) = lim u(t, s) , (9b) 

lim d s u(t, s) = lim d s u(t, s) , i G (0, 00) , (9c) 

and initial condition 

u(0,«)=uo(«)=uo(7(*)), sG(0,L). (9d) 
Hence, we obtain an arclength-parametric solution by Fourier series: 

00 1 r 

u(t,s)=e-» 2t cwe-^'e^)", ^ = - / u (s) e<^) ms ds . (10) 

m— — 00 q 

An arclength parametrization may not be given, so we also express (10) in terms of an arbitrary 
regular parametrization 7 : (— tt, tt) — > 5 \ {0}, f? — > 7(6'). The arclength parameter s = a{6) and 
the parameter 9 are related via the arclength function: 



a(9)= j W(r)\dr, (11) 



so that r y(a(9)) — 7(6*). As a matter of notation, we let u(t, 9) denote the parametric solution with 
respect to the parametrization 7(6), i.e. u(t, 9) = u(t, "/(9)). The relation between u and u is given 
by 

u(t, a(9)) = u(t, 7(a(0))) = u(t, 7 W) = &(*, 0) • 
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Consequently, by a change of variables we obtain the following series representation: 

oo 1 r 

u(t,8) = e""'* J] c m e~^ mH e l ^) ma{e) , c m = j- / u (6) e -*(*) mo W a '(0) dfl . (12) 

m— — oo _ 7r 

The parametrization that we use is 7(6*) = (1/2(1 + cos(0)), 1/4(1 + cos(0)) sin(0)) (derived in 
Section 4.1) and is not based on arclength, so we will work with (12). 



4 Regularization by Desingularization of the Domain 

4.1 Resolution of the Singularity in the Curve y 2 — x 3 — x 4 

In this section we resolve a singular point by a desingularization technique from algebraic geometry 
which is known as a blow-up (see for example, [17], [16], and [27]). For our curve S (see Figure 1) 
this will amount to viewing the cuspidal curve 5 as a two-dimensional projection of a non-singular 
curve S that lies in three-dimensional space. The blow-up process also appears in the catastrophe 
theory literature, see for example [2]. 

Recall, the cuspidal curve S is given by the zero level set of ip(x,y) :— y 2 — x 3 + x 4 which is 
singular at the origin. The Jacobian vanishes at the origin D(p(0,0) = (0,0), hence the equality 
Dip(0, 0) • v = holds for all «eK 2 . This means that the tangent space ToS is equal to all of M 2 
instead of being a one-dimensional subspace. 

Now, we will work with the tangent space TqS in an algebraic fashion by examining all straight 
lines passing through the origin. Specifically, we will consider 

y 2 - x 3 + x A = , (13) 

y--x = 0, (14) 
e 

where (14) describes a family of straight lines through the origin with different slopes z/e £ R. 
Figure 2 a) shows the situation for e = 1. 
Substituting (14) into (13), yields 



x 2 



{^)-x + x 2 ^ = 0, (15) 

y--x = 0. (16) 

We may apply Definition 1 to (15) by noting that the substitution step is equivalent to evaluating 
if on the straight line parametrized by (t,tz/e) 

.z\ .9 / / z \ - 2 



^M-j=^(-j -t + fj . (17) 

We see that ip(t, tz/e) as a polynomial in t has a root at t — with multiplicity 2 and this is true 
for all possible choices of z/e. This implies that any straight line passing through the origin is 
tangent to the singular curve S according to Definition 1. 

The key to the blow-up is to interpret the system of equations (15) and (16) in K 3 with z as 
the third variable. Each of the equations (15) and (16) defines a surface in R 3 and the intersection 
of the surfaces is a union of two curves. For a visual representation of this, see Figure 2 b). The 
first of these two curves is the z-axis, since if x = then (15) is already satisfied while (16) implies 
only y = 0, hence all triples (0, 0, t) solve the system. The second curve is obtained by considering 
the case where x 7^ and thus satisfies the system 

-x + x 2 = 0, (18) 

y--x = 0. (19) 
e 
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a) All straight lines through the origin are tangent to S (red) b) Intersection of the level-sets (13) and (14) (e = 1) 

z 




c) The original cuspidal curve (red) is the projection of the d) Intersection of the level-sets (18) and (19) (e = 1) 

desingularizcd curve (green) onto the surface z — 

Fig. 2: Visualization of the blow-up process for a parameter value of e = 1. a) shows the singular 
curve S (13) in red and some tangents (z is the slope) at the origin (14) in blue, b) visualizes the 
system of equations (13) (red surface) and (14) (blue surface) where z now represents the third 
dimension. The intersection of the two surfaces (green) contains the desingularized curve that we 
seek. Note that the intersection also contains the z-axis. c) shows the desingularized curve in green, 
detached from any surfaces, d) visualizes the system of equations (18) (red surface) and (19) (blue 
surface). In contrast to b), the red surface is now such that the intersection of the two surfaces 
(green) is exactly the desired desingularized curve. 



since we can divide by x. Note that the origin (0,0,0) is also a solution of (18) and (19) and 
thus belongs to both of these curves. Comparing (15) and (18) we see that the x 2 -factor which 
caused the singularity according to our discussion above is now gone. Hence, this second curve is 
the desingularization that we seek and we display it in Figure 2 c). 
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Finally, we complete the square in (18) with respect to x and obtain an equivalent desingular- 
ized system: 



& 



1 



ey — zx — . 



(20) 
(21) 



Figure 2 d) illustrates the system of equations (20) and (21) as the intersection of two surfaces. 
As can be seen from Figure 2 d) or the form of the equation, (20) defines an elliptical cylinder. 

The curve S e described by the system of equations (20) and (21) can be parametrized by 
7 e : [—ir, 7r) — > R 3 , where 





"1 0" 




%(0) := H e • 7(0) with H e := 


1 


, m ■= | 




e 





|(l + cos(0)) 
i(l + cos(0))sin(0) 
| sin(0) 



(22) 



The system of equations (20) and (21) defines a one-parameter family of desingularized curves 
S £ depending on e. The parameter value e = 1 corresponds to the desingularized reference curve 



S which is parametrized by 7 : [—it, tt) 



of (22). All further desingularizations S e of our 



one-parameter family are deformations of S where the linear deformation map is given in (22) by 
the matrix H e . 

As e > is reduced, the curve S e will become an increasingly accurate approximation of the 
original singular curve S. Indeed, So is a copy of S embedded in the ccy-plane of K 3 . In the situation 
where e = 0, the matrix Ho is the orthogonal projector onto the xy-plane, hence 70 parametrizes 
Sq and we obtain a parametrization 7 of the original singular curve S C M 2 by dropping the trivial 
z-component of 70. By reducing the domain of 7 to the open interval (— 7r,7r) we get the regular 
parametrization of the smooth part S \ {0} mentioned in Section 3.2. 

Finally, we note that the blow-up process provides us with a one-to-one correspondence between 
(x, y, z) S S e and (x, y) € <S*. The first map of this correspondence is clearly the projection 



S e ^ S , (x,y,z) -> (x,y) . 
The second map is the blow-up map, the inverse of this projection, and is given by 



S^S £ 



(0,0,0) 



x ^ 
x = 



(23) 



(24) 



4.2 The Regularized Problem 

The regularization is carried out by considering the same reaction-diffusion PDE on the new 
domain S s . Because S e is a smooth curve, we need not introduce an interior boundary. Instead, 
we have the following initial value problem (IVP) as the regularized problem 

d t v E = ^5 e w £ - ^ 2 v E , t € (0, 00) , (x, y,z) e S e , £ > (25a) 
v e (0,x,y, z) = u Q (x, y) , (x, y, z) € S e . (25b) 

The initial condition (25b) is constructed by lifting the initial condition of the original problem. 
This is formally done via the one-to-one correspondence between (x, y, z) £ S e and (x, y) G S. The 
projection map of (23) allows us to lift functions / defined on S to functions g defined on S e by 
/ -> .9, 9(x,y,z) := f(x,y). 

Because the solution v e of the regularized problem does not have the same domain as the 
solution u of the reference problem, we introduce an approximant u £ . The approximant u e will be 
defined on the original domain S and must tend to the desired solution u as e — > 0. Note that the 
blow-map of (24) induces a pull-down of functions g defined on S £ to functions / defined on S 
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by g — > /, f(x,y) := g(x,y,ey/x) (given that x ^ 0, otherwise we have /(0,0) = g(0,0,0)). We 
define the approximant u e as the pull-down of v e , i.e., 

w £ :(0,oo)xS\{0}^M, u £ {t, x, y) := v £ (t, x, y, e|) , (26) 

where i> e is the solution of (25a)— (25b). In Section 4.3 we will prove that for any time t > the 
family U e (t, ■ ) converges uniformly to u(t, ■ ) — the solution of the reference problem — as e tends 
to zero. 

The parametric solution of the regularized problem is easily constructed. The relevant parametriza- 
tion % of S £ is given by (22) and we denote the parametric solution by v e (t, 9) :— v e (t,%(9)). 
Employing the same method as in Section 3.2, we obtain 

v e (t,0)=e-»' t ]T c s , m e-(H) e <%) ma 'W , c e , m = ± / u (0) e<^) ma ^ a' e {6) d6 

m=-oc £ ^ 

— TV 

(27) 

where a £ (9) is the arclength function corresponding to %, i.e., 

o 

a £ {9):= J \%{r)\dr. (28) 

— 7T 

We will denote the length of the curve S s by L £ := a £ (ir). Finally, we note that setting the 
parameter value e to zero in (27) gives vo(t, 9) — u(t, 9), i.e. the parametric solution of the reference 
problem (12), even though the PDE (25a) is not defined for e = 0. 



4.3 Convergence Theory (e — > 0) 

We now turn our attention to proving the convergence of u s (t, ■ ) to u(t, ■). As a first step, we 
establish continuity of v £ in e = by the following lemma. 

Lemma 1 Let v £ (t, ■ ) denote the parametric solution of the regularized problem as given in (27) 
for parameter values e > 0. Let u(t, ■ ) denote the parametric solution of the reference problem as 
given in (12). Then for any fixed time t > 0, v £ (t, ■ ) converges uniformly to u(t, • ) as e — > 0. 

Proof We begin by estimating the difference a' e (9) — a' (9), where a £ (9) and Oq(9) = a{9) are the 
arclength functions from (28) and (11) respectively. 

< a' e {9) - a' (9) = a' e (9) - a' (9) = \%(9)\ - \%(9)\ < \%{9) - %{9)\ = < £ . (29 ) 

This estimate implies directly 

a' e {6) -a'{9) = 0(e) uniformly in 9 , 

o 

a e {9) - a{9) = J a' e (r) - a' (t) dr = 0(s) uniformly in 9 , (30) 

— TV 

L £ -L = a £ {n) - a(w) = 0(e) . 
From (30) we derive further implications which are dependent upon m but are uniform in 9: 
2itm 2nm 



L £ L 

2 



-(l + 0(e)) => 0> 



{-Lr)-{-^)= 0{mch (31) 



irma £ (9) 2-Kma{9) , „. ., 2ixma F (ff) 2 r Kma(ff) 

T = j^(l + 0(e)) => jf^- 1 ^=0{me) 
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The quantities defined below will be used later on 

2 /,-, \ 2 



a £: „j := exp 



(32) 

but we estimate them now. Using the results of (31) we see that a £ . m = 1 + (D(m 2 e) and that 
Se,m{9) = l + 0(me) for arbitrary but fixed m € Z as e — > 0. Hence, if M € N denotes an arbitrary 
but fixed bound on the indices to, we can say that 

a £ . m = 1 + 0{e) , for all \m\ <= M 
Se,m(0) = 1 + 0(e) , for all \m\ <= M , uniformly in 6. 

Furthermore, since 5 £jTn lies on the complex unit circle, we have 

%,m{6) - 1| < 2 , for all to G Z , uniformly in d. (34) 

We now proceed by examining the difference in the coefficients c £ rn given by (27) and c m given 
by (12): 

u (0) (e-<%) ma °W °M _ e -i(^)m-W ^ d0 . (35 ) 



•'S.m L-m 



Estimating the modulus |c e , m — c m |, we arrive at the following bound 

7T 

/ \Md)\ 



Cm I < 



— 7T 



°'<»> +0(£) )- W( ,)^ 



d0 



(36) 



< l|Mo||oo r l|a l|o ° / |l-<5 E , m (0)| d0 + O(e) 



Using (34) we observe that there is a bound B independent of to such that |c e m — c m \ < B. 
Moreover, the same bound holds for \c m \ and |c £iTri |, specifically \c m \ < B and |c SlTO | < B. Using 
the estimate (33) with our bound on m, we also find that 

|c E , m - c m \ = 0{e) , for all |m| <= M . (37) 



Next, we estimate the difference \v E (t, 9) — u(t, 8)\. Since t > is fixed the factors e - ^ 77 ^' m and 

e v^ 7 ' m decrease rapidly with increasing \m\. Hence for an arbitrary but fixed 5 > 0, we can 
choose M € N large enough so that the tails of \v s (t, 9)\ are small: 



-Af-l 

TT-) ma A8)\ 



m=M+l m=-oo 

oo -M-l oo r 

< £ Be -(!7) 2 »^ + £ S e -(^) 2 ™ 2 ' = 2i? ^ e^(^) 2m2t <^. 

m— Af+l m— — oo m— M+l 

A similar estimate holds for the tails of u(t,8). We note that the choice of M depends also on t. 
Because e virJ m is monotonically decreasing in to we can apply a standard integral test 

oo 00 / M \ 



m=M+l 



A/ 
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to obtain a suitable M. Now, using the estimates of the tails, we have: 

M 

\v e (t,6) - u(t,6)\ < K™ e-(^) 2 " l2t e ^m^) - Cm e -(*) 2 ™ 2 * e <¥)™*W\ + i 5 . 

m=—M 

(40) 

Rearranging the remaining summands, we arrive at the following inequality: 

M 2 2 

\v £ {t,9)-u(t,6)\ < J2 {\ce, m -c m \a e , m + \c m \\a e}m -6 £ , m {d)\} cOt) 2 ™ 2 * + -<5 . (41) 

m=-M 

Since the indices m are between — M and M, our bounds on a £ , m and S £ , m (9) as of (33) apply. 
This gives us the following estimate which is uniform with respect to 8 

\v e (t, 6) - u(t, 0)| < {0(e)(l + O(e)) + BO(e)} (2M + 1) + h < O(e) + 2 -5 . (42) 

Finally, we choose e small enough such that the O(e)-expression in (42) is lower than 5/3. This 
implies eventually 

||» e (i, •)-«(*, Oil- < 5 , ( 43 ) 
which proves the statement as S > was arbitrary. □ 

We now prove our main convergence result. 

Theorem 1 Let u e (t, •) be given by (26), i.e., as the pull-down of the solution v £ {t, ■) of the 
regularized problem (equations (25a)-(25b)J at time t. Let u(t, •) be the solution of the refer- 
ence problem (equations (8a)-(8d)J at time t. Then, for any fixed time t > 0, u £ {t, ■) converges 
uniformly to u(t, ■) as e — > 0. 

Proof We use the parametrizations 7 and 7 e from Section 4.1. By construction, the first two 
components are equal %j — 7;, i € {1, 2}, while the third component 7^3 is equal to 

7e.a - e- ■ (44) 
7i 

If we parametrize u £ (t, • ) by 7 then, using our convention u £ (t, 6) := u e (t, 7(6*)), we obtain 

u £ {t, •) =u e (t,j) = v £ ^,7i,72,e^ = v e (t,% A ,% j2 ,%,3) =v e (t,%) = v £ (t, •) . (45) 

Now, we appeal to Lemma 1 to obtain 

u £ (t,j)=u £ (t,-)~^u(t,-) = u(t,j) as £^-0, (46) 

where u(t, •) = u(t, 7) is the parametric solution of the reference problem. Since 7 is a regular 
parametrization of S \ {0} the statement follows from (46). □ 



4.4 Experimental Convergence Rate (e — > 0) 

Thus far we have shown that u £ (t, •) converges uniformly in space and pointwise in time to u(t, •) 
as e — > 0. For the purposes of developing a numerical method in later sections, we are interested 
in the rate at which u £ (t, ■) converges to u(t, ■). We investigate this by obtaining highly accurate 
approximations to a(0) and c m (as defined in (12)) as well as a £ {9) and c £ ^ m (as defined in (27)) by 
using Chebfun [ ] in MATLAB. We then employ the first equation of (12) and the first equation 
of (27) where we restrict the summation to — M < m < M with M chosen sufficiently large such 
that the bound e~ u M t \\u Q \\ oa is lower than machine accuracy. We then compute the Z°°-norm of 
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Fig. 3: Convergence results for various final times with various choices of e > 0. The initial 
condition is given by (47). 



the difference between u £ {t, •) and u(t, •). For our numerical experiments, we set fj, = 1 in (27) and 
(10), and choose the initial conditions: 

„ , „. exp(4cos 2 (f?)) . , , . 

6o(*, 0) = V 5Q V , O € [-7T, tt). (47) 

We perform a regression analysis on the linear part of the data; see Figure 3. Excluding the 
first six data points in the linear fit, we observe that the convergence rate of u e (t, •) to u(t, •) is 
E 1.85 m case _ However, we also note that in each case as the stop time tf decreases, the error 
curves in the double logarithmic plots of Figure 3 exhibit lower slopes for larger values of e. In 
addition, we ran similar tests with other periodic initial conditions possessing the same order of 
smoothness. In each case we found that the data exhibited similar convergence results with the 
convergence rate of u e (t, •) to u(t, •) ranging from e 183 to e 185 . 



4.5 Transformation of the Regularized Problem 

The regularized problem (equations (25a)-(25b)) is a one-parameter family of problems. In Sec- 
tion 5 these will be solved for various choices of e > as e —> by an embedding technique. 
However, as e — > 0, the domain S e approaches the singular domain S which makes it difficult to 
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solve the problem directly on the domains S s . We resolve this issue by transplanting the problem 
onto a fixed curve S = Si, in R 3 , and transforming the reaction-diffusion PDE to one involv- 
ing variable coefficients. To accomplish this task we introduce a new variable w £ via the linear 
diffcomorphism induced by the matrix H e of (22): 

w e (t,x) ■= v s (t,h e (x)) , h s :S^S e , h £ (x) = H e • x . (48) 

The pull-down given in (26) transforms accordingly to: 

u e (t, x, y) = w £ (t, x, y, |) , t e (0, oo) , (x,y)eS\ {0} . (49) 

Now, Lemma 2 below, tells us how to transform (25a) and (25b) into a problem involving w e . 

Lemma 2 Let r and r be two smooth curves embedded in M" which are diffeomorphic to each 
other under the linear diffeomorphism h : r — > r, h(x) = H • x. Let u : r — » K be a smooth 
function on r and let u : F — > K denote its pull-back u = u o h onto r . Then the pull-back of the 
Laplace- Beltrami operator Apuo h onto r transforms the following differential operator on F 

Aruoh= —divf,( ^Vf.-u) , (50) 

|H T| r V|H T| r J 

acting on u. In (50) T is (one of the two possible choices of) the unit tangent vector field of T. 

Proof According to Definitions 3 and 4, the Laplace-Beltrami operator on _T expands to 

A r u = div r (V r w) = trace(L»(PVu) P) , (51) 

where P(x) is the projector that projects onto the tangent space T^F . The right-hand side of (51) 
involves a C^-smooth local extension ue of u (cf. Definition 3), but we will simply use the identifier 
u for both the function and its extension. Next, we expand the right-hand side of (51) by using 
the product rule: 

A r u = trace (PD 2 uP + [D(Pw)] \ w=Va ■ P) . (52) 

Since r is a curve, the projector P(x) is simply the outer product P(x) = T(x) • T(x) T . We use 
the product rule again to expand the Jacobian of the vector Pw (where w independent of x) 

D(Pw) = D(TT T w) = DT{T T w) + Tw T DT . (53) 

With the last two results, the right-hand side of (52) becomes 

A r u = trace (D 2 uP) + trace (Z?T(T T Vw)P) + trace (T(V-u T L>T)P) 

= T T D 2 uT + (T T Vu)(T T DTT) + V/OTT = T T D 2 uT + \7u T DTT . ' ' ' ' 

In the last equality we have used that |T| 2 = 1 which implies T T Z?T = 0. Now, we pull back: 

A r u o h = T T o h (D 2 u) ohToh + (Vu T ) o h DT o h T o h . (55) 

The next step is to replace (Vw)o/i and (D 2 u)oh with expressions involving derivatives of u = uoh. 
We obtain these expressions from the chain rule. Moreover, the linearity of h, 

h(x) = H x , Dh = H , D 2 h = , (56) 

allows us to simplify in the following manner: 

Vu T = (V/) ohDh V/ o h = Vu T l-T 1 , (57) 

D 2 H = Dh T (D 2 u)ohDh D 2 uo h = hT T D 2 u H _1 . (58) 
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Finally, we must resolve Tok For our given h a regular parametrization 7 of r corresponds to a 
regular parametrization 7 of r via 7 = h(j). The linearity of h and the fact that the unit tangent 
vector field has parametric form T o 7 = 7'/|Y| yield the following transformation rule: 

y h 7' h y/iyi HT07 ht 

V; IYI |H 7 'I |H 7 'I/|7'I |HTo 7 | |HT| V ' 

Replacing the terms on the right-hand side of (55) with those of (57), (58), and (59), we arrive at 

T T L> 2 uf T ( f \ f 

A r u oh = = + Vu T D ^ . 60 

|H T| 2 \\H T\J |H T| 

The right-hand side of (60) involves only functions defined on r and can be simplified further to: 

Aruoh= — div^f i^Vpii I , (61) 

|H T| F V|H T| F J 

thus completing the proof. □ 

By Lemma 2, the regularized problem transforms to 

dtw £ = f3 E divg (p e V giu E ) ~ /i 2 w e , t € (0, 00), (x,y,z)<ES (62a) 
w e (0,x,y,z) =u (x,y) , (x,y,z)e§, (62b) 

where the non-constant coefficient function j3 e : S — ¥ M. is given by 

/8 B (x) = 1 , (63) 

|H £ -T(x)| 

and T : S —> M 3 denotes the tangent field of S. Notice that the new domain S is smooth and now 
independent of e. 



5 An Embedding Technique for the Numerical Solution 

In this section we solve the transformed regularized problem given by (62a) and (62b) without 
employing a parametrization; instead, we embed the problem itself into the ambient space M 3 . 
The transformation of Section 4.5 made the domain S independent of e. The advantage of this is 
that the choice of the computational grid is not affected by e, but we have traded this advantage 
for a nearly singular coefficient since /3 E (0) = 1/e. This makes the new problem stiff. 



5.1 Embedding and Discretization of the Problem 

Our numerical calculations will be carried out using the Closest Point Method [23,26,20], which 
is an embedding technique for solving PDEs posed on a given regular surface r of arbitrary 
codimension that is smoothly embedded in W 1 . We begin by reviewing some of the fundamental 
ideas behind the Closest Point Method. 
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5.1.1 Basics of the Closest Point Method 

The Closest Point Method represents a given smooth surface using a closest point function [21]. A 
closest point function cp : B(r) — > F maps every point x of a band B(r) (or tubular neighborhood) 
of r to a point cp(x) 6 J 1 and has the two properties: 

a) cp is a retraction, i.e., cp(x) = x if x £ r, 

b) cp is continuously differentiable with Z3cp(x) = P(x), i.e., the Jacobian of cp is equal to the 
projector P onto the tangent space of r whenever the point x is on the surface r. 

The properties can be used to establish the following two principles [21] which are fundamental 
to the Closest Point Method [26,20]: 

1. Gradient principle: Vpii(x) = V[u o cp](x), if x £ r, 

2. Divergence principle: divr <?(x) = div[g o cp](x), if x £ r, 

where u : r — >• R is a scalar surface function and g : r — > M. n is a surface vector field. The 
Gradient and Divergence principles allow us to replace surface intrinsic differentials with Cartesian 
differentials applied to the closest point extensions uocp and gocp. Since closest point extensions 
are defined on B(r) these principles when applied to a PDE posed on r yield an embedding 
equation posed on the band B(r). 

We now elaborate further on our construction of the closest point function. 

5.1.2 A Closest Point Function For the Desingularized Domain 

For the domain S of problem (62a)— (62b) we have an implicit description 

<p(x,y,z):=z 2 +(x-^J -~=0, (64) 

ip(x, y, z) := y — zx = (65) 

from Section 4.1. Following the approach of [21, Section 5], we construct — in a two-stage procedure — 
a closest point function cp : B(S) — > S directly using the given ip and (p. Let Sep and denote 
zero-level surfaces of tp and tp respectively. The intersection of these is exactly S. The first stage 
is to map a given point x 6 B(S) onto the surface S^. This is done by intersecting the trajectory 
£ of steepest decent, i.e., 

C' = -V^o£, £(0)=x (66) 
with the surface S^. This defines a unique point xGSf The second stage is then to map the new 
point x £ Sip onto S. Here, we intersect the trajectory r\ of 

7/ = -Vs^°r?, 7?(0) = x (67) 

with S. The second stage is well-defined because Vg^^ is always tangent to S^, hence the trajectory 
r] lies completely on the surface Sw, and evolves towards S. The second intersection is the point 
cp(x) £ S. This construction defines a differentiable retraction. The proof that this defines a 
closest point function can be found in [21]. 

Now we discuss how to define the band B(S) which forms the domain of cp. Since VV> does not 
vanish we can solve for £ for arbitrary initial points x. However, V<p vanishes on the straight line 
{(1/2, y,0) : y £ M}. The intersection of this line with is the point (1/2,0,0). Consequently, 
Vs^c/? will vanish only at (1/2,0,0) and thus we cannot solve for rj by (67) if x = (1/2,0,0). By 
ODE (66) we observe that the set of points that are mapped to x = (1/2, 0, 0) in the first step is 
given by the trajectory £(r) = (cosh(r)/2,r, — sinh(r)/2), rel. This implies that the maximal 
band is B max (S) = R 3 \ £(R). In practice we define and use a smaller band that is contained in 
B ma x(S)] for example, the following band is also admissible under the above criteria: 

B(S) := | (x, y, z) £ R 3 : ^5 tp(x, y, z) 2 + ^(x, y, z)* < \ J . (68) 

This choice of band is depicted in Figure 4. 
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Fig. 4: A possible choice of the band B(S) (red) defined according to (68). The blue solid line is 
the curve S, while the black dashed line is £, the set of points non-projectable by the closest point 
function cp as constructed in Section 5.1.2. The maximal possible band is B max (S) = R 3 \ £(R). 



5.1.3 The Embedding Equation 

The embedding equation is obtained by applying the Divergence and Gradient principles to the 
original problem. For this purpose we define the linear extension operator E by Eu(x.) := uocp(x). 
The benefit of this is that the Divergence and Gradient principles may be written as Ediv^g = 
EdivEg and EV§Eu, see [22]. We now apply the extension operator to (62a) and (62b) and 
arrive at the following embedded evolution equation: 

d t Ew £ = ELEw £ - ^Ewe , t € (0, oo) , x = (x, y, z) e B(S) , 

(69) 

Ew E (0, x) = w (cp(x)i, cp(x) 2 ) , x = (x, y, z) G B(S) , 
where the operator L is defined by 

Lf := & divE (&V/) , &(x) := & o cp(x). (70) 

Finally, we write the embedded evolution equation (69) in the equivalent form 

w E (0, x) = u (cp(x)i, cp(x) 2 ) (71a) 
d t w £ — {EL — [i 2 I)w e , te(0,oo),x=(x,y,z)eB(S), (71b) 
w E = Ew e . (71c) 

Here, we have replaced Ew e with w £ and used the fact that E is idempotent [21,22,12]. The 
operator / in (71b) denotes the identity operator. 

5.1.4 Discretization of the Embedding Equation 

We discretize the embedded problem given by (71a)-(71c) in a similar manner to [26,22], i.e. (71b) 
is discretized in time with either a forward Euler (Algorithm 1) or a backward Euler (Algorithm 2) 
time step (of step size r) followed by an extension in order to satisfy (71c). For the discretization in 
space, we use a uniform Cartesian grid Gh of mesh width h on the box [—0.5, 1.5] x [—1, 1] x [—1, 1]. 
Computations are performed within the banded grid Gh D B(S) and we denote a point in the grid 
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Algorithm 1 Explicit Iteration 

™h,e '■= u o(cp(xh)i,cp(x h ) 2 ) 
for n = 0,1,2, ... do 

v>h ■= (Ih + T{E h L h - fJ. 2 I h )) ■ w 

end for 



Algorithm 2 Implicit Iteration 

a h,e : = u o(cp(x h )i,cp(x /l ) 2 ) 
for n = 0,1,2, ... do 

Wh ■= {h ~ r(E h L h - A* 2 /ft))" 1 ■ t&Jk 

Kt 1 ■= E h ■ ™h 

end for 



Table 1: Explicit and implicit closest point iteration for problem (71a)-(71c). 



by Xft. We further use the subscript h for spatially discrete functions (vectors) and operators 
(matrices). The superscript n is used for temporally discrete functions and refers to time nr 
where r is our time step-size. The spatially discrete operators in Algorithms 1 and 2 are Ih, Eh, 
Lh- Here, Ih is simply the identity matrix while Eh is an interpolation matrix. Because cp(x^) 
is in general not a grid point, we approximate the function value of a closest point extension 
Ew(xh) —wo cp(xfc)- This is done by performing a tri-cubic interpolation with the discrete 
function data Wh on the 4 x 4 x 4-point neighborhood surrounding the point cp(x^) [26,20]. 
Finally, we expand the operator L from (70) in terms of partial derivatives and obtain Lh as a 
composition of matrices: 

Eh ■= B h ■ (D Xjh E h B h D Xth + Dy jh E h B h D Vth + D z hE h BhD z h ) , (72) 

where Bh is a diagonal matrix representing the factor /3 6 . The matrices E x h, D y h, D z h ap- 
proximate the partial derivatives d x , d v , d z by central differences, and Eh is again the extension 
matrix. 



5.1.5 Consistency and Stability 

In [22], it is shown that the truncation error of both the explicit and the implicit scheme is of order 
0(h 2 + h 4 /t + t). Thus we choose the time step-size to be r = 0(h 2 ) to obtain good accuracy. This 
gives us an overall truncation error of 0(h 2 ) and, assuming stability, an anticipated convergence 
rate of 0(h 2 ). The following heuristics guide our choice of time step-size. Notice that problem 
(62a) is stiff: the function is largest at the origin (3 £ (0) = 1/e, e <C 1, hence the spatial operator 
behaves like 

(3 £ dWg (f3 e W giv e ) — /J 2 w e s» —^Agiv e — /i 2 ui £ as x > . (73) 

Since G cxp i := Eh(I + r(Lh — ^ 2 I)), the iteration matrix of the explicit iteration (Algorithm 1), 
is obtained from the forward Euler time-stepping scheme, we arrive at the following time-step 
restriction: 

t = C cxp i • e 2 • h 2 . (74) 

As we are interested in small parameter values e<1, this implies a very large number of time 
steps to reach a given stop time. In contrast to this, the implicit iteration with Gi mp i := Eh(I — 
r(EhLh — /i 2 /)) -1 , which is based on backward Euler, is unconditionally stable (regarding e) and 
can therefore be based on accuracy concerns alone, i.e., 

r = C lmpl • h 2 . (75) 

In our experiments, the time step-size r is chosen according to (74) with C eX pi = 1/4 in the explicit 
case and according to (75) with C; mp i = 1 in the implicit case. 
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Explicit iteration, Tstop = 0.001 Implicit iteration, Tstop = 0.001 




10" 3 10" 2 10"' 10"' 10 2 10"' 

h h 



Fig. 5: Double logarithmic plot of the errors in Algorithm 1 (left column) and in Algorithm 2 
(right column) measured in the Z°°-norm at stop times t = 0.001 (first row) and t = 0.1 (second 
row). The initial condition is given in (76) and the reaction parameter is \i = 1. 



5.2 Numerical Examples 

We run Algorithms 1 and 2 for different choices of e £ {0.5,0.05,0.005} and measure errors by 
comparing the results to the independently obtained parametric solutions of Section 4.4. The 
initial function is 

exp (4 ■ (2.x - l) 2 ) 

M x >y) = — - — ^ -> ( 76 ) 

which corresponds to (47). The parameter fi in the reaction term is set equal to /i = 1. 

Figure 5 shows the errors as a function of the mesh width h in double logarithmic plots. The 
left column gives results for the explicit iteration while the right column gives results for the 
implicit iteration. The first row gives the error at a stop-time of tf — 0.001, and the second row 
gives the error at a stop-time of tf = 0.1. In all cases, the error was measured in the ?°°-norm over 
the parameter domain 9 € [— tt, tt). Examining the right column, we observe a convergence rate 
of 0(h 2 ), and that the error constant is relatively independent of the parameter e which caused 
the stiffness. For the explicit iteration, we observe an C(/i 2 )-convergence for fixed e, however the 
error constant grows as e decreases. The number of time steps needed for the explicit iteration to 
reach a certain stop time becomes excessively large for small values of e (compare to (74)). This 
explains why the red curve is missing in the bottom left plot of Figure 5. Clearly, the implicit 
iteration is much better suited for our stiff problem since the time step-size is independent of e. 

Summarizing, we see that the reference problem (equations equations (8a)-(8d)) can be solved 
by combining the blow-up method with an embedding technique that can handle arbitrary co- 
dimensions such as the Closest Point Method. 
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6 Extensions and Generalizations 

6.1 Regularization by Desingularization of the Domain 

In Section 4 we regularized the reaction-diffusion problem posed on the cuspidal curve shown in 
Figure 1. The main step was the desingularization of the PDE domain via the resolution of the 
cusp singularity. In general we consider reaction-diffusion PDEs posed on real algebraic surfaces 
with finitely many isolated singularities and want to apply the same regularization technique. 

In order to resolve a single singularity of a real algebraic surface embedded in E™ we blow up 
E™ at the point p £ E" of singularity. As in Section 4.1 we investigate the tangent space T p S of S 
at p by adjoining equations that describe all straight lines through p. Any straight line through p 
can be parametrized by p + tz with a non-zero vector z £ E™ and parameter t £ E. If x is a point 
on this straight line then x — p and z must be linearly dependent. Consequently, the matrix 



xi - pi Zi 



Pn %n 



(77) 



has rank < 1 and all its 2 x 2-minors vanish. The vanishing 2 x 2-minors yield an implicit description 
of the straight line through p in direction z by 

(xi - Pi)zj - (xj - p 3 )zi = , 1 <= i < j <= n . (78) 

In (78) there is plenty of redundancy since there are n(n— 1)/2 equations in x to describe a straight 
line in R™. Hence, we can choose a subset of n — 1 linearly independent equations which we adjoin 
to our original system of polynomial equations that define the surface S. 

The so-called blow-up of E ra at the point p £ E" is given by viewing (78) as a system of 
equations in variables x and z. As on the one hand z is non-zero and on the other the length of 
z is irrelevant, it is standard in algebraic geometry to consider z as an element of the projective 
space P™ _1 , i.e., (x, z) £ E" x P" _1 . For more details on blow-ups the reader is referred to [17,16, 
27]. Since we want to blow up in E" x E n_1 , we use one component of z as the parameter e, i.e., 
z — (e,z) £ E x E n_1 , and view (78) as a system of equations in variables x and z. This yields a 
one-parameter family of blow-ups in E™ x E n . 

The final steps to desingularize are: substitute (78) into the original system of polynomial 
equations, pull out the factors that caused the singularity and neglect them. In the end we obtain 
a system of polynomial equations which describes a family S £ , e > 0, of smooth surfaces that 
correspond to the original surface S by projection, i.e., e = 0. Given several singularities, we 
resolve one at a time by the same procedure. We will demonstrate this procedure later with 
more examples. In addition, we point out that blow-up calculations can be automatized with the 
computer algebra system SINGULAR [ ]. 

This procedure for the resolution of singularities generalizes that of Section 4.1. Now, we would 
like to regularize, as in Section 4.2, by considering the same reaction-diffusion problem on S £ . But 
this is not possible if the blow-up in E™ x E" _1 does not preserve the topology or connectivity of 
the original domain of the PDE. We now give two examples of what can go wrong. 

The first example is the cardioid (see left plot of Figure 6) which is given by 

(x 2 +,/) 2 -x(x 2 + y 2 )-^=0. (79) 
The procedure above or that of Section 4.1 yields the following desingularization: 



z 2 \ 1\ 2 1 / z 2 



z 

y x = 



(80) 
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Fig. 6: Left: cardioid in K 2 given by (79). Right: cardioid embedded in xy-plane of R 3 (red) and 
its desingularization (x > green, x < blue) given by (80) with e = 1/4. The blow-up procedure 
tears the curve apart. 




Fig. 7: Left: figure-eight curve in R 2 given by y 2 — x 2 + x A = 0. Right: figure-eight curve embedded 
in xj/-plane of M 3 (red) and its desingularization (green) given by (82) with e = 1/4. The blow-up 
procedure pulls the origin apart. 



From this description we derive the following parametrization j e : [— n, ir) — > W 3 , y e = ((1 + 
cos(0)) cos(0)/2, (1 + cos(0)) sin(0)/2, etan(0)) of S £ . Here, the blow-up in M 2 x R 1 tears the curve 
apart, since | tan(0)| tends to infinity as approaches ±7r/2 (see right plot of Figure 6). The blow- 
up changes the topology compared to the original curve which was a closed loop. Consequently, 
we cannot use the same reaction-diffusion problem on S £ to regularize the problem. 

As second example consider the figure-eight curve (see left plot of Figure 7) given by 

y 2 - x 2 + x 4 = . (81) 

The desingularization yields 




A parametrization -% : [— ir,ir) —> M 3 is given 7 6 = (sin(0), cos(6>) sin(0), e cos(6>)). It is clear from 
this parametrization that the z-axis intersects with the regularized curve S e twice. Indeed, the 
blow-up procedure regularized the curve by pulling it apart at the origin as shown in the right 
plot of Figure 7. 

This situation is very different from the cardioid because here it depends on the boundary 
conditions if we obtain a consistent regularization or not. The original reaction-diffusion problem 
consists of a PDE on the figure-eight curve excluding the origin, together with boundary conditions 
at the origin. If we impose boundary conditions such that the origin is a cross-junction then 
resolution of the singularity changes the topology since the green curve in the right plot of Figure 7 
exhibits a different connectivity. In this case considering the same reaction-diffusion PDE on S e will 
not give a consistent regularization of the problem. If we are given boundary conditions such that 
the connectivity matches that of the desingularization, then the same reaction-diffusion problem 
on S £ gives us the desired regularization of the problem. 
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u @ t=0 u @ t=0.001 




Fig. 8: Reaction-diffusion on the e = 0.5-blow-up of the singular surface S given by y 2 +z 2 = x 3 — x A . 
The images show the pull-down u £ of the solution v £ of the IVP (84) onto the singular surface 
S. The left image shows the initial state u £ at time while the right image shows u £ at time 
t f = 0.001. 



A sufficient condition for a consistent regularization is that the blow-up procedure does not 
change the topology of the surface. We conjecture — since there is evidence in our examples — 
that this condition is satisfied if for every singular point p there is a hypcrplane passing through p 
which docs not intersect the surface at a different point, i.e., the surface is contained in a half-space 
defined by this hyperplane through p. 



6.2 Surfaces with an Isolated Singularity 

So far, we have considered curves with a cusp singularity, however, we are also interested in 
applying our strategy to surfaces with similar singularities. Consider the singular surface S given 
by y 2 + z 2 = x 3 — x (see Figure 8). This is the surface of revolution derived from our familiar curve 
(6). In order to desingularize S we follow the same steps as in Section 4.1. The only difference 
is that a description of all straight lines that pass through the origin requires now two equations 
since S is embedded in IR 3 . Similar to the case of our cuspidal curve, we obtain a one-parameter 
family of desingularized surfaces S £ which are embedded in M 5 . These are given by the system 

c n (83) 
ez — i]x — . 

We remark that all surfaces S s are linear deformations of S = Si, and that the deformation matrix 
is H e = diag{l, 1, l,e,e}. 

The numerical results displayed in Figure 8 are for the reaction-diffusion problem 

d t v e = A§ e v s - v s , t € (0,oo) , x e S £ , u e (0,x) = u (x) , (84) 

where S £ is the regular two dimensional surface embedded in K 5 given by (83) with e = 1/2 and 
x := (x, y, z, £, rf). The initial condition is 

, , exp (2 • cos(4arccos(2a; — l)) 2 ) , v 

«o(x) = ? 5 }±± ■ (85) 
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Fig. 9: The red curve is given by y 2 = x 3 (l — x) 3 and possesses two isolated singularities. 

To obtain a domain amenable to numerical computation we transform, as in Section 4.5, via 
io E (x) = v £ (H £ x) to an equivalent problem on S. This yields a new domain independent of e and 
the variable coefficient PDE problem: 

d t w £ — trace (Cj D§ (C E V gW E )) — w E , t£(0,oo), x€S, w s (0, x) = uo( x ) • (86) 

This transformation is a generalization of Lemma 2 to surfaces of arbitrary dimension. The oper- 
ators Dg and Vg denote the Jacobian and the gradient intrinsic to S and the coefficient matrix 
C e is given by the expression 

C e = H e Pg (Hj 2 - H- 2 N (iV T H i T 2 iV)- 1 ^ T H7 2 ) . (87) 

The matrix TV is the transpose of the Jacobian of system (83), and its columns span the normal 
space. 

We solved the regularized problem (86) using the explicit closest point iteration (Algorithm 1). 
Our computations were performed on a banded grid around the surface Scl 5 . The banded grid 
is a subset of the virtual 21 x 21 x 21 x 21 x 21 Cartesian grid defined on the box [—0.3, 1.2] x 
[-0.75,0.75] x [-0.75,0.75] x [-0.75,0.75] x [-0.75,0.75] and contains 0.05 % of the points of full 
virtual grid. Figure 8 shows u e at times t = and t = 0.001, where u e denotes the pull-down of 
v e onto the singular surface S: 

u e (t,x,y,z) = v s (t,x,y,z,e-,s-) = w e (t,x,y,z, -) . (88) 
V xx) \ xx/ 

This higher dimensional numerical result illustrates that the ideas presented are applicable to 
certain algebraic surfaces of arbitrary dimension and codimension where the surface has isolated 
singular points. Of course there are limitations regarding the type of singularity, for example, 
a singular point must be such that its blow-up is a smooth surface in M™ which has the same 
topology as the original singular surface (cf. Section 6.1). 

This example also illustrates that the desingularized surface is embedded in a high dimensional 
space in general. This side effect leads us to pay close attention to numerical efficiency. For example, 
a computational band should be used that is optimal or nearly optimal; see e.g. [20, Appendix A]. 

Recall that for the problem of the cuspidal curve we proved in Section 4.3 that the pull-down 
u £ converges to the solution u on the original cuspidal curve without the singular point. Such 
a convergence investigation for surfaces would be an interesting subject for future work. Note, 
however, that a different analysis may be required due to the inherently more complex geometry 
of two dimensional surfaces over one dimensional curves. Assuming convergence of u e given by 
(88) as e ~ > 0, we would guess that the limit u is a solution of the problem 

d t u — A s u - u , t £ (0,oo) , x = (x,y, z) e S\ {0} , u(0,x) = tio(x) , (89) 

on the domain S \ {0} where S = {(x, y, z) £ E 3 : y 2 + z 2 = x 3 — x 4 }. 
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6.3 Curves and Surfaces with Multiple Singularities 



When finitely many isolated singularities arise, our procedure is to resolve one singularity at a 
time. We illustrate this generalization by resolving the two isolated singularities of the curve shown 
in Figure 9 which is given by 

y 2 - x 3 (l - xf = . (90) 

The singular points are the origin (x,y) = (0,0) and the point (x, y) — (1,0). Starting out with 
the point (x, y) — (1, 0), we adjoin an equation describing the straight lines passing through (1, 0). 
This yields 

y 2 - a; 3 (l - a;) 3 = , 
V [x - 1) = . 

El 

The next steps agree exactly with Section 4.1. Specifically, we replace y in the first equation of 
(91) with (x — l)z/e, we factor out (1 — x) 2 , and we desingularize by neglecting this factor. This 
yields 



2 

z ' 



El 



- x 3 {l -x) = 0, 

(92) 

y--(x-l) = 0. 

El 



In the new system (92), only the first equation is singular. Indeed, it is exactly our example from 
Section 4.1. We therefore adjoin an equation for the straight lines passing through the origin 
(x,z) — (0,0) as in Section 4.1. This yields the desingularized system 



1\ 2 1 



x — 



eie 2 ; v 2/ 4 ' 

y--(x-l)=Q, (93) 

El 

w n 
z x = . 

E2 

This system defines a curve in four dimensional space, (x, y, z, w) G K 4 , which can be parametrized 
by 7(0) = ((1 + cos(£?))/2,-sin(6>) 3 /8,£isin(6l)(l + cos(0))/4, £ X £ 2 sin(0)/2). The first two com- 
ponents of 7 give a parametrization of the original singular curve. Formally, we can introduce a new 
parameter Si for each singularity. But looking at the deformation matrix H £lje2 = diag{l, 1, Ex, £i£2} 
we observe that it is sufficient to keep only parameter e% = e and set the other equal to £2 = 1- 
The deformation matrix H e = diag{l, l,e, e} will still give us the projection in the case of £ = 0. 
Accordingly, the desingularized system is 



(?) 



2 / lV 1 



e) ' \ X 2) 4 



z ( x -D = o, (94) 



z — wx = 



Finally, we remark that the dimension of the embedding space increases only by one per 
resolved singularity. After the resolution of the first singularity we obtained system (92) in M 3 , so 
one might have expected that for the resolution of the second singularity we had to adjoin two 
equations for the description of a straight line in IR 3 (cf. Section 6.2). That would have yielded 
a system of four equations in five variables. Interestingly, the resolution of a singularity does not 
change the number of variables in the first equation. Thus, every single resolution operation can 
be performed by considering only the first equation in two dimensional space. 
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Analogously, for a two dimensional algebraic surface in K 3 with several isolated singularities we 
adjoin two equations each time we resolve a singularity. For example, consider the surface obtained 
from revolving the curve shown in Figure 9 around the £-axis: 

y 2 + z 2 - x 3 (l - x) 3 = . (95) 

The resolution of the singularities (one singularity at a time) gives the desingularized system 



b\ 2 ( 1\ 2 1 



'a 

X ~ 2j " 4 ' 
y- a (x-l) = 0, 

z - -(x- 1) = , 
e 

a — cx = , 
b - dx = . 

This is a system of five equations involving seven variables (x, y, z, a, b, c, d) £ R 7 as opposed to a 
system of seven equations in nine variables. 



7 Conclusion 

This paper presents a framework and numerical embedding method for evolving reaction-diffusion 
equations on algebraic surfaces with isolated singularities. We explored and proved several analyt- 
ical results in both geometry and geometric PDEs in Section 4. In particular, we used the blow-up 
technique from algebraic geometry to resolve singularities in a domain, proved the convergence 
of the solutions for the reaction-diffusion equation on S e to the solution of the reaction-diffusion 
equation on S, and constructed a variable-coefficient PDE on S that corresponds to a PDE posed 
on the domain S e . Together these results provide a solid mathematical foundation from which we 
built a numerical technique. 

In constructing our numerical method, we first applied Lemma 2. This lemma allows us to pose 
a variable-coefficient PDE on a smooth and numerically well-resolved curve S that corresponds 
to the reaction-diffusion PDE on S e . We then applied the Closest Point Method to evolve the 
variable-coefficient reaction-diffusion equation. Our numerical results are stable and demonstrate 
that the method obtains second order accuracy in space. 

In Section 6 we discussed extensions as well as limitations of our approach. In particular, we 
saw that preservation of the topology plays an important role. Notably, the convergence theory 
of Section 4.3 applies to all closed curves with finitely many singularities as long as the blow-up 
procedure preserves the topology of the curve. This is because the preservation of the topology 
implies that (29) in the proof of Lemma 1 holds (with a different constant). 

As desingularization of the PDE-domain proved to be useful for regularization we are in- 
terested in combining embedding techniques such as the Closest Point Method with desingular- 
ization other than blow-up. This might help to reduce some of the limitations of our method 
that have been identified in this paper. For example, the problem of topological change with 
the cardioid as discussed in Section 6.1 is resolved if we use the curve parametrized by j e — 
((1 + cos(0)) cos(0)/2, (1 + cos(0)) sin(0)/2, e sin(0)). This curve is suitable for regularizing the 
cardioid since it corresponds to the cardioid by projection, is smooth, and has the same topol- 
ogy as the cardioid. Besides further extending our study of surfaces, the investigation of other 
desingularization techniques will be an interesting subject for future work. 
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